Neophobia analysis

Setup

Load and install R packages:

packages <- c('lmerTest', 'lme4', 'ggplot2', 'tidyverse', 'readxl', 'purrr', 'performance', 'emmeans', 'MASS')

# Check if each package is installed; if not, install it
for (pkg in packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg)
  }
  library(pkg, character.only = TRUE)
}

Data Loading

Grab neophobia_data.csv from processed_data directory:

script_dir <- dirname(rstudioapi::getActiveDocumentContext()$path)
# Set the working directory to the script's directory
setwd(script_dir)

data <- read.csv("processed_data/neophobia_data.csv", row.names = 1)

Quick Data Inspection

Summary statistics of the data:

# Quick look at structure
str(data)

'data.frame':   536 obs. of  15 variables:
 $ Bird_ID         : chr  "BB_BB" "BB_BB" "BB_BB" "BB_BB" ...
 $ Trial           : int  1 2 3 4 5 6 7 8 1 2 ...
 $ Enclosure       : chr  "B8" "B8" "B8" "B8" ...
 $ Object          : chr  "novel" "novel" "control" "control" ...
 $ Context         : chr  "group" "individual" "group" "individual" ...
 $ Latency_to_enter: num  2.2 6.23 2.27 1.57 1.07 ...
 $ Latency_to_Eat  : num  2.9 13.83 2.2 3.73 2.73 ...
 $ Zoi_duration    : num  183.7 91.9 85.4 203.3 371.5 ...
 $ Object_Type     : int  2 4 3 3 1 5 3 3 5 1 ...
 $ GroupID         : chr  "7A" NA "7A" NA ...
 $ NestID          : chr  "2024VM08" "2024VM08" "2024VM08" "2024VM08" ...
 $ group_dummy     : int  1 0 1 0 1 0 1 0 1 0 ...
 $ ind_dummy       : int  0 1 0 1 0 1 0 1 0 1 ...
 $ Object_contrast : num  0.5 0.5 -0.5 -0.5 0.5 0.5 -0.5 -0.5 0.5 0.5 ...
 $ Context_contrast: num  0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 ...

summary(data)

   Bird_ID              Trial       Enclosure            Object            Context          Latency_to_enter  Latency_to_Eat   
 Length:536         Min.   :1.00   Length:536         Length:536         Length:536         Min.   :  0.800   Min.   :  1.792  
 Class :character   1st Qu.:2.75   Class :character   Class :character   Class :character   1st Qu.:  1.467   1st Qu.:  2.800  
 Mode  :character   Median :4.50   Mode  :character   Mode  :character   Mode  :character   Median :  1.880   Median :  3.633  
                    Mean   :4.50                                                            Mean   :  8.380   Mean   : 42.457  
                    3rd Qu.:6.25                                                            3rd Qu.:  2.600   3rd Qu.:  5.684  
                    Max.   :8.00                                                            Max.   :600.000   Max.   :600.000  
  Zoi_duration     Object_Type      GroupID             NestID           group_dummy    ind_dummy   Object_contrast Context_contrast
 Min.   :  0.00   Min.   :1.000   Length:536         Length:536         Min.   :0.0   Min.   :0.0   Min.   :-0.5    Min.   :-0.5    
 1st Qu.: 38.93   1st Qu.:2.000   Class :character   Class :character   1st Qu.:0.0   1st Qu.:0.0   1st Qu.:-0.5    1st Qu.:-0.5    
 Median : 76.35   Median :3.000   Mode  :character   Mode  :character   Median :0.5   Median :0.5   Median : 0.0    Median : 0.0    
 Mean   :128.20   Mean   :2.787                                         Mean   :0.5   Mean   :0.5   Mean   : 0.0    Mean   : 0.0    
 3rd Qu.:162.03   3rd Qu.:4.000                                         3rd Qu.:1.0   3rd Qu.:1.0   3rd Qu.: 0.5    3rd Qu.: 0.5    
 Max.   :598.13   Max.   :5.000                                         Max.   :1.0   Max.   :1.0   Max.   : 0.5    Max.   : 0.5    

Data Preprocessing

Missing Group ID Values

The GroupID values are currently coded as NA during individual trials. Assign the most frequent group for each bird:

most_frequent_group <- function(group_ids) {
  group_ids <- group_ids[!is.na(group_ids)]
  if (length(group_ids) == 0) return(NA)
  names(sort(table(group_ids), decreasing = TRUE))[1]
}

most_frequent_group_per_bird <- data %>%
  group_by(Bird_ID) %>%
  summarize(most_frequent_group = most_frequent_group(GroupID)) %>%
  ungroup()

data <- data %>%
  left_join(most_frequent_group_per_bird, by = "Bird_ID") %>%
  mutate(GroupID = ifelse(is.na(GroupID), most_frequent_group, GroupID)) %>%
  dplyr::select(-most_frequent_group)

Adjusting Trial Numbers, so the baseline is the first trial rather than 0.

Adjust the Trial variable by subtracting 1:

data$Trial <- data$Trial - 1

Statistical Modeling

1) Latency to Enter

The full model as described in the RR, includes Object_contrast, Context_contrast, and Trial, and a complex random effect structure:

enter_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial + 
                (1 | NestID) + 
                (-1 + group_dummy | GroupID) + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data)

boundary (singular) fit: see help('isSingular')

summary(enter_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +  
    (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy +      group_dummy | Bird_ID)
   Data: data

REML criterion at convergence: 5712.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4529 -0.1408 -0.0422  0.0686 10.2384 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 Bird_ID  ind_dummy   1.198e+03 34.6124      
          group_dummy 2.852e-02  0.1689  1.00
 NestID   (Intercept) 0.000e+00  0.0000      
 GroupID  group_dummy 0.000e+00  0.0000      
 Residual             2.268e+03 47.6258      
Number of obs: 536, groups:  Bird_ID, 67; NestID, 43; GroupID, 15

Fixed effects:
                                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                       17.4867     4.3340 279.0124   4.035 7.06e-05 ***
Object_contrast                    5.1073     4.1143 464.9955   1.241  0.21510    
Context_contrast                 -12.0511     5.8852  79.6165  -2.048  0.04389 *  
Trial                             -2.6020     0.9052 471.1684  -2.874  0.00423 ** 
Object_contrast:Context_contrast  -9.7670     8.2292 464.9964  -1.187  0.23588    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_ Trial 
Objct_cntrs  0.005                     
Cntxt_cntrs -0.356  0.000              
Trial       -0.731 -0.007  0.007       
Objct_cn:C_ -0.010  0.000  0.000  0.013
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

check_model(enter_model)

The random effect structure seems to be too complex for the amount of data. However, the variance for the (Intercept) under NestID is 0 and the variance for group_dummy under GroupID is 0. This suggests that both nest as differences between groups contribute minimally to the variance in the outcome model. Let’s drop both effects.

enter_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data)

boundary (singular) fit: see help('isSingular')

summary(enter_model2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +      (-1 + ind_dummy + group_dummy | Bird_ID)
   Data: data

REML criterion at convergence: 5712.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4529 -0.1408 -0.0422  0.0686 10.2384 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 Bird_ID  ind_dummy   1198.0356 34.6127      
          group_dummy    0.0286  0.1691  1.00
 Residual             2268.2114 47.6257      
Number of obs: 536, groups:  Bird_ID, 67

Fixed effects:
                                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                       17.4867     4.3340 279.0054   4.035 7.06e-05 ***
Object_contrast                    5.1073     4.1143 464.9960   1.241  0.21510    
Context_contrast                 -12.0511     5.8852  79.6162  -2.048  0.04389 *  
Trial                             -2.6020     0.9052 471.1689  -2.874  0.00423 ** 
Object_contrast:Context_contrast  -9.7670     8.2292 464.9969  -1.187  0.23588    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_ Trial 
Objct_cntrs  0.005                     
Cntxt_cntrs -0.356  0.000              
Trial       -0.731 -0.007  0.007       
Objct_cn:C_ -0.010  0.000  0.000  0.013
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

check_model(enter_model2)

Seems like the interaction between object and context is non-significant, let’s drop it. Given the non-normal distribution, let’s logtransform the data as well:

enter_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data)

boundary (singular) fit: see help('isSingular')

summary(enter_model3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Latency_to_enter) ~ Object_contrast + Context_contrast +      Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
   Data: data

REML criterion at convergence: 1123.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9335 -0.4596 -0.0811  0.3538  6.9301 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Bird_ID  ind_dummy   0.307211 0.55427      
          group_dummy 0.008857 0.09411  1.00
 Residual             0.386408 0.62162      
Number of obs: 536, groups:  Bird_ID, 67

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)        1.29304    0.06319 188.76075  20.464   <2e-16 ***
Object_contrast    0.10177    0.05370 465.99922   1.895   0.0587 .  
Context_contrast  -0.19502    0.07774  76.68931  -2.508   0.0142 *  
Trial             -0.13771    0.01179 468.78462 -11.679   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_
Objct_cntrs  0.004              
Cntxt_cntrs -0.458  0.000       
Trial       -0.653 -0.007  0.007
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

check_model(enter_model3)

There are still issues with the data distribution, let’s remove birds that did not interact and boxcox_transform:

data_enter <- data %>% filter(Latency_to_enter != 600)

boxcox_transform <- boxcox(lm(Latency_to_enter ~ Object_contrast + Context_contrast + Trial, data = data_enter))

best_lambda <- boxcox_transform$x[which.max(boxcox_transform$y)]

data_enter$Latency_to_enter_trans <- (data_enter$Latency_to_enter^best_lambda - 1) / best_lambda

enter_model4_boxcox <- lmer(Latency_to_enter_trans ~ Object_contrast + 
                              Context_contrast + Trial + 
                            (-1 + ind_dummy + group_dummy | Bird_ID), 
                            data = data_enter)
summary(enter_model4_boxcox)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter_trans ~ Object_contrast + Context_contrast +      Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
   Data: data_enter

REML criterion at convergence: -115.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.00860 -0.62234  0.05268  0.62220  2.42212 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Bird_ID  ind_dummy   0.023257 0.15250      
          group_dummy 0.007021 0.08379  0.84
 Residual             0.037006 0.19237      
Number of obs: 533, groups:  Bird_ID, 67

Fixed effects:
                   Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)        0.686462   0.020720 165.959664  33.130   <2e-16 ***
Object_contrast    0.007148   0.016674 395.856100   0.429    0.668    
Context_contrast  -0.028455   0.020248  64.608541  -1.405    0.165    
Trial             -0.051896   0.003669 424.210940 -14.146   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_
Objct_cntrs  0.007              
Cntxt_cntrs -0.301 -0.002       
Trial       -0.623 -0.009  0.016

check_model(enter_model4_boxcox)

Comparing the different models, the boxcox transfromed comes out as best:

enter_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial + 
                (1 | NestID) + 
                (-1 + group_dummy | GroupID) + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data_enter)

boundary (singular) fit: see help('isSingular')

enter_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data_enter)

boundary (singular) fit: see help('isSingular')

enter_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data_enter)

boundary (singular) fit: see help('isSingular')

enter_model4_boxcox <- lmer(Latency_to_enter_trans ~ Object_contrast + 
                              Context_contrast + Trial + 
                            (-1 + ind_dummy + group_dummy | Bird_ID), 
                            data = data_enter)


anova(enter_model, enter_model2, enter_model3, enter_model4_boxcox)

refitting model(s) with ML (instead of REML)

Data: data_enter
Models:
enter_model3: log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
enter_model4_boxcox: Latency_to_enter_trans ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
enter_model2: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
enter_model: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
                    npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)
enter_model3           8  959.8  994.0  -471.91    943.8                     
enter_model4_boxcox    8 -127.5  -93.2    71.73   -143.5 1087.3  0           
enter_model2           9 5148.8 5187.3 -2565.41   5130.8    0.0  1          1
enter_model           11 5152.8 5199.9 -2565.41   5130.8    0.0  2          1

2) Latency to Eat Model

A similar model is built for Latency_to_Eat, again incorporating interaction terms and random effects:

eat_model <- lmer(log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial + 
                    (1 | NestID) + 
                   (- 1 + group_dummy | GroupID) + 
                   (- 1 + ind_dummy + group_dummy | Bird_ID),  
                 data = data)

boundary (singular) fit: see help('isSingular')

summary(eat_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +  
    (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy +      group_dummy | Bird_ID)
   Data: data

REML criterion at convergence: 1636.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1780 -0.5149 -0.1108  0.1927  4.6862 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Bird_ID  ind_dummy   0.78398  0.8854       
          group_dummy 0.08989  0.2998   1.00
 NestID   (Intercept) 0.00000  0.0000       
 GroupID  group_dummy 0.00000  0.0000       
 Residual             1.00669  1.0033       
Number of obs: 536, groups:  Bird_ID, 67; NestID, 43; GroupID, 15

Fixed effects:
                                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                        1.89202    0.10740 165.73052  17.616  < 2e-16 ***
Object_contrast                    0.74392    0.08668 465.00491   8.583  < 2e-16 ***
Context_contrast                  -1.00808    0.11239  90.13817  -8.969 3.94e-14 ***
Trial                             -0.04274    0.01899 466.49134  -2.251   0.0248 *  
Object_contrast:Context_contrast  -0.85108    0.17337 465.00516  -4.909 1.27e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_ Trial 
Objct_cntrs  0.004                     
Cntxt_cntrs -0.434  0.000              
Trial       -0.619 -0.007  0.008       
Objct_cn:C_ -0.008  0.000  0.000  0.013
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

check_model(eat_model)

Similarly, the variance for the (Intercept) under NestID is 0 and the variance for group_dummy under GroupID is 0. Let’s drop those:

eat_model2 <- lmer(log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial + 
                   (- 1 + ind_dummy + group_dummy | Bird_ID),  
                 data = data)
summary(eat_model2)
check_model(eat_model2)

There are still some issues with the normality. Let’s remove non participating birds, and try the same boxcox tranformation:

boxcox model seems to fit data best

eat_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial + 
                (1 | NestID) + 
                (-1 + group_dummy | GroupID) + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data_eat)

boundary (singular) fit: see help('isSingular')

eat_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data_eat)

boundary (singular) fit: see help('isSingular')

eat_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data_eat)

boundary (singular) fit: see help('isSingular')

eat_model4_boxcox <- lmer(Latency_to_eat_trans ~ Object_contrast + 
                              Context_contrast + Trial + 
                            (-1 + ind_dummy + group_dummy | Bird_ID), 
                            data = data_eat)


anova(eat_model, eat_model2, eat_model3, eat_model4_boxcox)

refitting model(s) with ML (instead of REML)

Data: data_eat
Models:
eat_model3: log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
eat_model4_boxcox: Latency_to_eat_trans ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
eat_model2: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
eat_model: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
                  npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)
eat_model3           8  690.1  723.9  -337.03    674.1                     
eat_model4_boxcox    8 -794.8 -761.0   405.39   -810.8 1484.8  0           
eat_model2           9 4137.0 4175.0 -2059.50   4119.0    0.0  1          1
eat_model           11 4141.0 4187.5 -2059.50   4119.0    0.0  2          1

3) Time spent in the ZOI

The full model as described in the RR, includes Object_contrast, Context_contrast, and Trial, and a complex random effect structure:

summary(enter_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +  
    (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy +      group_dummy | Bird_ID)
   Data: data_enter

REML criterion at convergence: 4110

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.7468 -0.0958 -0.0198  0.0611 19.4975 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 Bird_ID  ind_dummy   1.196e+02 10.93520     
          group_dummy 3.693e-04  0.01922 1.00
 NestID   (Intercept) 0.000e+00  0.00000     
 GroupID  group_dummy 0.000e+00  0.00000     
 Residual             1.722e+02 13.12072     
Number of obs: 506, groups:  Bird_ID, 67; NestID, 43; GroupID, 15

Fixed effects:
                                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                        5.2873     1.2911 247.2473   4.095 5.72e-05 ***
Object_contrast                   -0.8899     1.1804 435.9808  -0.754   0.4513    
Context_contrast                  -1.7996     1.7837  70.9317  -1.009   0.3164    
Trial                             -0.5942     0.2595 438.0996  -2.289   0.0225 *  
Object_contrast:Context_contrast   1.9819     2.3611 435.9820   0.839   0.4017    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_ Trial 
Objct_cntrs  0.057                     
Cntxt_cntrs -0.433 -0.040              
Trial       -0.722 -0.038  0.031       
Objct_cn:C_ -0.057 -0.071  0.045  0.041
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

The random effect structure seems to be too complex for the amount of data. However, the variance for the (Intercept) under NestID is 0 and the variance for group_dummy under GroupID is 0. This suggests that both nest as differences between groups contribute minimally to the variance in the outcome model. Let’s drop both effects.

zoi_model2 <- lmer(Zoi_duration ~ Object_contrast * Context_contrast + Trial + 
                (-1 +  ind_dummy + group_dummy | Bird_ID), 
              data = data)
summary(enter_model2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +      (-1 + ind_dummy + group_dummy | Bird_ID)
   Data: data_enter

REML criterion at convergence: 4110

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.7472 -0.0958 -0.0198  0.0611 19.4973 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 Bird_ID  ind_dummy   1.196e+02 10.93614     
          group_dummy 3.683e-04  0.01919 1.00
 Residual             1.722e+02 13.12060     
Number of obs: 506, groups:  Bird_ID, 67

Fixed effects:
                                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                        5.2873     1.2912 247.2062   4.095 5.72e-05 ***
Object_contrast                   -0.8899     1.1804 435.9864  -0.754   0.4513    
Context_contrast                  -1.7996     1.7838  70.9217  -1.009   0.3165    
Trial                             -0.5942     0.2595 438.1047  -2.289   0.0225 *  
Object_contrast:Context_contrast   1.9818     2.3611 435.9877   0.839   0.4017    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_ Trial 
Objct_cntrs  0.057                     
Cntxt_cntrs -0.433 -0.040              
Trial       -0.722 -0.038  0.031       
Objct_cn:C_ -0.057 -0.071  0.045  0.041
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

check_model(enter_model2)

Seems like the interaction between object and context is non-significant, let’s drop it. Given the non-normal distribution, let’s logtransform the data as well:

summary(enter_model3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Latency_to_enter) ~ Object_contrast + Context_contrast +      Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
   Data: data

REML criterion at convergence: 696

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4872 -0.5604 -0.0780  0.4192  8.9279 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Bird_ID  ind_dummy   0.083162 0.28838      
          group_dummy 0.009698 0.09848  1.00
 Residual             0.196081 0.44281      
Number of obs: 506, groups:  Bird_ID, 67

Fixed effects:
                   Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)        1.143564   0.043897 223.101604  26.051   <2e-16 ***
Object_contrast    0.006807   0.039692 438.913282   0.171    0.864    
Context_contrast  -0.074761   0.046032  93.240157  -1.624    0.108    
Trial             -0.114456   0.008698 438.594691 -13.159   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_
Objct_cntrs  0.052              
Cntxt_cntrs -0.324 -0.045       
Trial       -0.711 -0.035  0.038
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

There are still issues with the data distribution, let’s boxcox_transform:

Comparing the different models, the logtransfromed comes out as best:

anova(zoi_model, zoi_model2, zoi_model3, zoi_model4_boxcox)

refitting model(s) with ML (instead of REML)

Data: data
Models:
zoi_model3: log(Zoi_duration) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
zoi_model2: Zoi_duration ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
zoi_model4_boxcox: Latency_to_zoi_trans ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
zoi_model: Zoi_duration ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
                  npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)
zoi_model3           8 1342.4 1376.2  -663.21   1326.4                     
zoi_model2           9 6314.2 6352.3 -3148.11   6296.2    0.0  1          1
zoi_model4_boxcox    9 1416.4 1454.4  -699.18   1398.4 4897.9  0           
zoi_model           11 6313.3 6359.8 -3145.63   6291.3    0.0  2          1

So we end up with these models:

# Time near object
summary(zoi_model3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Zoi_duration) ~ Object_contrast + Context_contrast + Trial +      (-1 + ind_dummy + group_dummy | Bird_ID)
   Data: data

REML criterion at convergence: 1343.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.95388 -0.63743 -0.01154  0.62983  2.78860 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Bird_ID  ind_dummy   0.2019   0.4493       
          group_dummy 0.2805   0.5296   1.00
 Residual             0.6831   0.8265       
Number of obs: 506, groups:  Bird_ID, 67

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)        4.38246    0.09114 174.24034  48.086   <2e-16 ***
Object_contrast    0.06530    0.07394 438.97006   0.883    0.378    
Context_contrast   0.01085    0.07456 357.90017   0.145    0.884    
Trial              0.02323    0.01618 438.49089   1.436    0.152    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Objct_ Cntxt_
Objct_cntrs  0.045              
Cntxt_cntrs  0.036 -0.047       
Trial       -0.636 -0.035  0.042
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Combined Latency Model

This model considers the combined effect of latency and different contrasts:

latency_model <- lmer(log(Latency) ~ eat_vs_leave_contrast * Object_contrast * Context_contrast + Trial + 
                    (- 1 + group_dummy | GroupID) + 
                    (- 1 + ind_dummy + group_dummy + eat_vs_leave_contrast | Bird_ID), 
                  data = data)

Estimated Marginal Means (EMMs)

Compute and compare estimated marginal means for the eat_model:

emm_eat <- emmeans(eat_model, ~ Object_contrast * Context_contrast)
pairs(emm_eat, adjust = "bonferroni")

Zone of Interest (ZOI) Model

Finally, a model is built for the Zoi_duration variable:

zoi_model <- lmer(Zoi_duration ~ Object_contrast * Context_contrast + Trial + 
                    (1 | NestID) + 
                    (1 + group_dummy | GroupID) + 
                    (1 + ind_dummy + group_dummy | Bird_ID), 
                  data = data)
summary(zoi_model)
check_model(zoi_model)

Additional Analysis

As mentioned in the research report (RR), we may also consider analyzing latencies combined in a multivariate model:

multivariate_model

(Note: Details of the multivariate model should be added based on the specific analysis required.)

Conclusion

This document provides a comprehensive workflow for analyzing neophobia data, from data loading and preprocessing to statistical modeling and result interpretation. Further steps can be added as needed to refine and extend this analysis.